library(Amelia)
library(Zelig)

dataset<-data.frame(read.csv("dem_demise.csv"))
dataset.new<-dataset[c("country","year","gdp_annual_growth","gdp_growth_5_year","gdp_baseline","log_gdp_per_capita","log_ppp_per_capita","log_gdp_baseline","gov_expend",
"log_inflation","trade_percent","remittance","nat_resource","net_bilateral",
"net_multilateral","inequality","imr","ethnic","independence","pres","liberal",
"caseid","demage","demduration","reversed","xconst",
"region","pre1980","lat_am","e_eur","ss_africa","asia","postcol")]

# run Amelia to impute missing data
a.out<-amelia(x=dataset.new,m=7,ts="year",cs="country", ords=c("xconst","nat_resource","postcol","liberal"),
idvars=c("reversed","caseid","pre1980","ethnic","independence","lat_am","e_eur","ss_africa","asia"),interecs=TRUE)
model_columns <- c("demduration", "reversed", "caseid", "gdp_growth_5_year", "log_gdp_baseline", "log_ppp_per_capita", "log_inflation", "gov_expend", "trade_percent", "nat_resource", "ethnic", "imr", "inequality", "independence", "postcol", "liberal", "net_bilateral", "net_multilateral", "lat_am", "e_eur", "ss_africa", "asia", "xconst", "pres")

# get out only the columns we want from each imputation
slice <- a.out$imputations
for (x in seq(1,7)) {
	slice[[x]] <- slice[[x]][model_columns]
}

# latin american observations
subset.lat_am <- slice
for (x in seq(1,7)) {
	subset.lat_am[[x]] <- subset.lat_am[[x]][subset.lat_am[[x]]["lat_am"] == 1,]
	subset.lat_am[[x]] <- subset.lat_am[[x]][-which(apply(subset.lat_am[[x]], 1, function(y) {(is.na(y))})),]
	subset.lat_am[[x]] <- na.omit(subset.lat_am[[x]])
}

# eastern european observations
subset.e_eur <- slice
for (x in seq(1,7)) {
	subset.e_eur[[x]] <- subset.e_eur[[x]][subset.e_eur[[x]]["e_eur"] == 1,]
	subset.e_eur[[x]] <- subset.e_eur[[x]][-which(apply(subset.e_eur[[x]], 1, function(y) {(is.na(y))})),]
	subset.e_eur[[x]] <- na.omit(subset.e_eur[[x]])
}

# sub-saharan african observations
subset.ss_africa <- slice
for (x in seq(1,7)) {
	subset.ss_africa[[x]] <- subset.ss_africa[[x]][subset.ss_africa[[x]]["ss_africa"] == 1,]
	subset.ss_africa[[x]] <- subset.ss_africa[[x]][-which(apply(subset.ss_africa[[x]], 1, function(y) {(is.na(y))})),]
	subset.ss_africa[[x]] <- na.omit(subset.ss_africa[[x]])
}

# asian observations
subset.asia <- slice
for (x in seq(1,7)) {
	subset.asia[[x]] <- subset.asia[[x]][subset.asia[[x]]["asia"] == 1,]
	subset.asia[[x]] <- subset.asia[[x]][-which(apply(subset.asia[[x]], 1, function(y) {(is.na(y))})),]
	subset.asia[[x]] <- na.omit(subset.asia[[x]])
}

# hazard ratio a / b
hazard_ratio <- function(alpha,betas,a,b) {
	(exp(b %*% betas)^alpha) / (exp(a %*% betas)^alpha)
}

# draw our betas from zelig's model fit
draw_hazard <- function(z.out, baseline, increased) {
	beta.draws <- mvrnorm(2000, mu=c(coefficients(z.out), log(z.out$scale)), Sigma=z.out$var)

	p.ests <- c()
	alpha_draws <- c()
	beta_draws <- rbind()
	for (i in 1:2000) {
		thisDraw <- beta.draws[i,]
		beta <- thisDraw[1:(length(thisDraw)-1)]
		shape <- thisDraw[length(thisDraw)]
		p.ests[i] <- hazard_ratio(1/exp(shape), beta, increased, baseline)
		alpha_draws <- c(alpha_draws, 1/exp(shape))
		beta_draws <- rbind(beta_draws, beta)
	}
	list(p.ests, alpha_draws, beta_draws)
}

for (current_sub in list(subset.lat_am, subset.ss_africa, subset.asia)) {

# run zelig to get the weibull model
z.out <- zelig(Surv(demduration, reversed) ~ gdp_growth_5_year + log_gdp_baseline + log_ppp_per_capita + log_inflation + gov_expend + trade_percent + nat_resource + ethnic + imr + inequality + independence + xconst * pres + postcol + liberal + net_bilateral + net_multilateral, model="weibull", robust=TRUE, cluster="caseid", data=current_sub)

dim(z.out) <- length(z.out)

# find all the medians
medians <- apply(current_sub[[1]], 2, function(x) {median(x)})

# set the high
xhigh.out <- setx(z.out)
xhigh.out <- xhigh.out[-19]
for (x in colnames(xhigh.out)[-1]) {xhigh.out[x] <- as.numeric(medians[x])}
xhigh.out["xconst:pres"] <- as.numeric(xhigh.out["xconst"]) * as.numeric(xhigh.out["pres"])
#xhigh.out["log_ppp_per_capita"] <- log(1501)
#xhigh.out["ethnic"] <- 0.75
#xhigh.out["pres"] <- 0
#xhigh.out["xconst"] <- 7
#print(summary(whatif(data=current_sub[[1]][whatif_columns], cfact=xhigh.out[-1][-18])))

# set the low
xlow.out <- setx(z.out)
xlow.out <- xlow.out[-19]
for (x in colnames(xlow.out)[-1]) {xlow.out[x] <- as.numeric(medians[x])}
xlow.out["xconst:pres"] <- as.numeric(xlow.out["xconst"]) * as.numeric(xlow.out["pres"])
#xlow.out["log_ppp_per_capita"] <- log(501)
#xlow.out["ethnic"] <- 0.25
#xlow.out["pres"] = 0
#xlow.out["xconst"] <- 3
xlow.out["net_bilateral"] = 0

vals <- c()
alpha_vals <- c()
beta_vals <- rbind()
for (i in seq(1,7)) {
	results <- draw_hazard(z.out[[i]], as.matrix(xlow.out), as.matrix(xhigh.out))
	vals <- c(vals, results[[1]])
	alpha_vals <- c(alpha_vals, results[[2]])
	beta_vals <- rbind(beta_vals, results[[3]])
}

vals <- sort(vals)
alpha_vals <- sort(alpha_vals)

# mean and stddev of alpha (shape) parameter
print(mean(alpha_vals)); print(sqrt(mean(alpha_vals^2) - mean(alpha_vals)^2));

# mean and stddev of betas
print(apply(beta_vals, 2, mean)); print(sqrt(apply(beta_vals^2, 2, mean) - apply(beta_vals, 2, mean)^2));

# confidence intervals of hazard ratio
print(vals[700]); print(median(vals)); print(vals[13300])
}







